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Numerical data on the probability distribution of the equilibrium relaxation time of the 
Sherrington-Kirkpatrick model are obtained by means of dynamical Monte Carlo simulation, for 
several values of the system size A*' and temperature T. Proper care is taken that the thermal fluc- 
tuations on the relaxation time estimates are totally negligible compared to the disorder induced 
fluctuations. The probability distribution of In r — In r scales with the scaling variable A'^^^^ (Tc — T) 
strengthening the behef that Inr oc N^^^ in the whole spin glass phase. 
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The equilibrium dynamics of the Sherrington-Kirkpatrick (SK) model remains a subject of much interest. The 
^ . standard picture of the spin glass phase of this model is that of a complex hierarchical free energy landscape, with 
/~s ' many valleys that correspond to pure or metastable states. In the thermodynamic limit, both the number of valleys 
, and the height of the typical free energy barrier between two valleys go to infinity. Accordingly the relaxation time 
of the system diverges as the number of spins N goes to infinity. 



The behavior of the equilibrium relaxation time of the model t with the system size N has been studied by analytical 
methods p'-S] , direct Monte Carlo simulation of the Metropolis dynamics of the model 041] , and indirect determina- 
. tion of the largest barrier height [Tj] . There are reasonable indications that below Tc the disorder averaged relaxation 
I time behaves according to Inr B/T oc as TV — ?> oo, with B the largest barrier height and an exponent V' — 1/3, 
^ for both binary and Gaussian distributions of the couplings. The behavior of the relaxation time of this model with 
^ the system size N has also been studied in the aging (non equilibrium) regime ^ . 



A different numerical approach to this problem has been used recently by Monthus and Garel [Ifl]- (This method 
has already been used in |TT| for the 2d Ising model and in Q for the SK model.) They use the well know mapping 
C ' (see e.g. chapter 4 of [l^) of the master equation for the Monte Carlo dynamics onto an Schrodinger equation in 
configuration space with some quantum Hamiltonian Hj (where J stands for a particular disorder configuration). The 
(~| ground state of "Hj has zero energy, corresponding to the equilibrium stationary state of the Monte Carlo dynamics. 
O . The next eigenvalue is the inverse of the largest relaxation time tj of the Monte Carlo dynamics of the original model 
, (the so called exponential relaxation time of the dynamics [13 )■ The problem of determining tj is thus reduced to the 

' problem of finding the lowest eigenvalues of a real symmetric sparse matrix (of size 2^ x 2^), which can be obtained 
fT^ ■ with high accuracy using a standard computer routine. The process has to be repeated for a large number of disorder 
^ ' configurations J. 

QQ ■ Compared to the direct Monte Carlo method, the method of Monthus and Garel has two clear advantages: i) It is 
, not affected by thermal noise; ii) the long tail of the probability distribution P/v(lnr) can be easily sampled. Indeed, 
^3 provided a good starting point is guessed, the convergence of the eigenvalue search is very fast. The method is however 
limited to very small system sizes, indeed the analysis of [lo| relies on systems with 6 < A'^ < 20, to be compared 
with the nine years old direct simulations of [6|] where 64 < A < 1024. This is in principle a strong limitation for a 
I model that is critical in the whole low temperature phase with slowly decaying power law finite size corrections. For 
T-H , example for such small sizes the shape of the probability distribution of the order parameter P{q) is strongly affected 
^ ' by finite size effects and is quite different from the textbook shape of the infinite volume limit. 

The results of Monthus and Garel for the SK model with Gaussian couplings at temperature T ~ Tc/2 = 0.5 can be 
summarized as follows: The disorder averaged logarithm of the largest relaxation time behaves according to 



Inr -B/T oc A^''' V 0.33 , (1) 

N^oc 



on small systems ( A^ < 20) already. The probability density function of In r scales like 



/, N 1 p./lnr — lnr\ 
PA^(lnr) = ^p(^^), (2) 



using the measured values of = (Inr)^ — (Inr)'^ and Inr, with an A" independent P(.). Monthus and Garel are 
not able to measure the width exponent i/jwidth defined by A cx TV'/'midtfi (-which means that it is crucial to use the 
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measured values of A and Inr in equation [2l at least for small systems), but make a tentative indirect estimate 
'4'width — 0.26 from their measurement of the tail exponent rj (defined later in equation [5]) and an assumption about 
the disorder configurations that dominate the tail of P{x) for x » 1. In [3] already, the quoted value of ijjwidth ~ 0.25 
is lower than 1/3. We remark however that the figure 1 of [?[ shows systematic errors as big as 0.05 in the value of 
-0 from fits of data with 128 < N < 1024, which are blamed on finite size effects. The agreement between the results 
of [l3| and for tpwidth (that is harder to measure than is accordingly somewhat surprising. 

The purpose of this note is to investigate the distribution P/v(lnT) using the direct Monte Carlo simulation method, 
which gives access to much larger system sizes than the method of (lOj . and without the strong assumptions about 
the dynamics of the multi canonical algorithm made in We take 1024 disorder samples (with binary couplings) 
and measure for each sample the dynamical overlap 



'Zj(i) = ^^'7«(io)fT*(t + to) , (3) 



averaging over along a long trajectory starting from a well equilibrated spin configuration |18j . i.e. we measure 
the thermal averaged < qj{t) > at equilibrium. In practice a chain of lO'* Metropolis sweeps (with sequential site 
update Is*!) was generated for two independent copies of the system (two clones), starting from two independent 
spin configurations, for each disorder sample. Measurements were made every 4 sweeps (namely the average was 
done over the values to = 0,4,8,...). It would be a waste of CPU time to measure the value of qj{t) for every 
(integer) values of t. It was measured for the following values: [1-20] with lag 1, [22-40] with lag 2, [44-80] with 
lag 4, .... Altogether qj{t) was measured for 194 values of t, up to a maximum value of t^ax = 3670016. For 
N — 512, a smaller chain of 4 10^ Metropolis sweeps was generated, and qj{t) was measured for 183 values of t up to 
a maximum value of tmax = 1703936 only. The relaxation time tj is defined by the condition qj{Tj) = \/ < q^ > j/2, 
where the mean value of the static overlap squared < q^ >j has been measured at equilibrium in the same disorder 
sample. Note that the ratio qj{t)/ y/< q^ >j is dimensionless and is accordingly [l^ a function of t/rj. Namely 
qj{t)l \/< > J = Fj{t/T), with some Fj{.) that is a continuously decreasing function of its argument. The precise 
analytical form of Fj{.) is irrelevant. We have checked that, disorder sample by disorder sample, the difference 
between qj{t) measured using clone one and qj{t) measured using clone two is so small that both give in most cases 
indistinguishable results for tj. (The worst case is for = 64 and T = 1 with a relative error of 0.07 for a particular 
bad disorder sample, namely the observed discrepancy is always smaller than 0.07, when not exactly zero, for more 
interesting values of T and N.) This shows that the plus or minus one standard deviation estimates of qj{t) give 
the same estimate for tj to an excellent precision, and that accordingly the thermal noise is negligible. That the 
thermal errors are negligible can alternatively be seen as follows: At fixed J, with P independent measurements of 
In Tj, namely Inr^^ In Tj'^-', . . ., one has the elementary unbiased estimator of Sj, the thermal statistical error on Inr, 

given by the expression Sj = 1/{P — ^)(^^/ PJ2i'=ii^''^'''j^)'^ ~ i^/^J^iLi • our simulation the two clones 

provide two fully independent measurements of In tj, and we have accordingly the unbiased (but noisy) estimator 
Sj — l/4(lnTy'' — InTj'^'')^. As mentioned before, the ratio Sj/hiTj is very small for all disorder sample J , system 
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size A^ and temperature T. We note that the disorder averaged y Sj is less than one percent of the median of Inr for 
all values of A^ and T (For obvious reasons, we consider only the values of A^ and T such that the median of r is less 
than tmax, and the average of Sj is restricted to the disorder samples for which r < tmax)- The disorder average of 
is even smaller. 

The fact that the thermal noise is negligible is an essential condition in order to obtain meaningful estimates for 
the probability distribution of the relaxation time. The various parameters (the total run length, the window of 
measurement tmax, and the lag between two successive measurements) have been chosen empirically, they are such 
that the thermal noise is negligible (as we just said) , the CPU times spent in Monte Carlo updates and measurements 
are balanced and the program fits inside the computer memory. No attempt was made to optimize theses parameters, 
it should in principle be done separately for each disorder sample, each size N and each temperature T . We remark 
that the time scales considered in Q are different from the one considered here, indeed considers time scales (called 
Ti, T2 and Ts in [^) that can be defined from the time decay of qj{t) measured with a single starting point on the 
one hand, and time scales (called and Tg2 in Q) defined from the time decay of the disorder averaged < q{t) > (or 
< q{ty >) on the other hand . The time scales Ti, T2 and have both thermal and disorder fluctuations, whereas 
Tq and rg2have no disorder fluctuations by construction. The time scales considered in the present note are disorder 
dependent with negligible thermal noise contamination. The purpose of [6] was indeed to show that all time scales in 
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FIG. 1: (Color on line) Scaling of the probability distribution P]v(lnr) at fixed temperature T — 0.6. 



the SK model behave according to Inr oc N^^^. 

We have data for N = 64, 128, 256 and 512, with temperatures between T — 0.4 and T = 1.1 with steps of 0.1 (the 
critical temperature of the model is — 1). Figure [1] shows a scaling plot of PN{lnT)N^^^ as a function of Inr/iV^/'^ 
for T = 0.6. The relative statistical error on the value of Pn{-) inside a bin is 1/^/Q, with Q the number of data 
points inside the bin. Figure [1] show good scaling, confirming that Inr has a behavior compatible with ^ = 1/3, 
namely (In this formula P{.) is an N independent function, related to the function P(.) introduced above). 



extending the scaling of the probability distribution found in from iV < 20 to TV < 512, with both Inr and 
A explicitly proportional to N^^^. Data taken at other values of the temperature (below Tc) show similar scaling. 
Depending on the temperature our sampling of the tail of the probability distribution is limited by the number of 
Monte Carlo sweeps performed (this is the case at low temperature) or by the number of disorder samples used (this 
is the case at higher temperature). Since we are plotting the logarithm of the (histogrammed) probability distribution 
of the logarithm of the relaxation time, enlarging substantially the data range in figure [T] would require a huge increase 
of the computational effort. 

Implicit in figure [T] is that the width A scales with the same exponent 1/3 as the mean. This can be checked directly: 
Following 3 we have computed the median M{N) of the distribution of Inr and a width W{N) defined arbitrarily 

Iw(N) PN{^n.T) dlnr = 0.30. The use of the median of the distribution instead of the average, and of a width 
defined from the quantiles of the distribution has the advantage that the later is only defined if rj < tmax for all 
disorder samples whereas the former requires that at least one half of the disorder samples satisfy this bound. The 
computational gain is enormous with a distribution that has a very long tail towards large values of r. The ratio 
W{N)/M{N) as a function of N is shown in figure [5] for several temperatures (without statistical errors). It shows a 
weak N dependence, consistent with the same scaling for the width and the disorder average of the distribution. (As 
a function of T however the ration W{N)/M{N) decreases slightly as T grows.) Indeed a plot of Pjv(lnr)W^(A^) as a 
function of (Inr — M{N))/W{N) shows the same scaling as figure [TJ 

Our conclusion for the behavior of the distribution of Inr with N is that, going to quite larger system sizes confirm 
the results of Monthus and Garel. With larger systems however, there is no need anymore to use the measured mean 
value and width of the distribution, in scaling plots. Assuming that both scale like iV^/^ just does it. Comparing 
to previous numerical determinations of the exponent t/;, the fact that the whole distribution, and not only the 
median, scales with the exponent = 1/3 strengthens the conclusion that the exponent is indeed ip — 1/3. 

The behavior of the function P{x) for large values of x defines the tail exponent ry, 

\nP{x) cx -x" , (5) 
that is 7] ~ 1.36 according to Monthus and Garel [lOj. The value 77 = 1 would imply a linear slope in figure [1] whereas 
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FIG. 2: (Color on line) The ratio of the width W{N) of the distribution divided by the median M{N), as a function of A'^ for 
several temperatures. 
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FIG. 3: (Color on line) Scaling at fixed size A'' of the probability distribution P(lnr — Inr). Here A'^ — 256. 



the value 77 — 1.36 would imply a slight downwards curvature. Both values are clearly compatible with our data, that 
unfortunately do not sample deep enough inside the tail of P/v(lnT) to allow a meaningful estimate of 77. 

Since we have data for several values of the temperature, we can see if our probability distributions scale with the 
temperature also. We have seen that, at fixed temperature, the distribution of the largest barrier scales like iV^/'^. 
Since the Sherrington-Kirkpatrick model is a mean field model, the scaling combination is N^^^'^'^^p^Tc — T), where 
dup = 6 is the upper critical dimension of the theory and v = 1/2 [l^, |T^. The scaling of the probability distribution 
of the relaxation time has been known for some time for the spherical Sherrington-Kirkpatrick, and is of this form 
indeed [i], 



7V1/3 

lnrj~^(r,-T)R, (6) 

with R a random variable (independent of N and T). Here N"^^^ R is the difference between the two largest 
eigenvalues of the coupling matrix of the model. Mathematical proof of this scaling behavior and the expression for 
the probability distribution of R can be found in ■ 

In figure [3] we check the hypothesis that indeed In rj = Inr + (Tc — T)N^/^Ii/T with R a random variable (with a 
distribution that does not depend on N or T). This figure shows reasonable scaling, even if the temperatures are not 
so close to Tc- Indeed the scaling quality deteriorates if one adds points closer to Tc- This seems paradoxical, the 
likely explanation is that tj is not a pure exponential and that sub-leading power law contributions to tj becomes 
important close to Tc where the dominant term in the exponential goes to zero, since at Tc one has t oc N^^'^^p. In 
order to obtain a good scaling at fixed TV it is crucial to consider the distribution of In tj — In r as is done in figure |31 
There is no such need at fixed T like in figure [T] The likely explanation is that Inr behaves like N-^^^ but not like 
(Tc — T)/T, or that sub-leading power law contributions to tj are important. 
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Figure m shows the ratio A = (lnr)2/lnr , where the overhne is here the usual arithmetic average (using the median 
would give a trivial result). It should be a constant, independent of N and T, if tj was proportional to (Tc — T)N^^^. 
Small scaling violations are visible, the data for increasing number of spins N seem to converge toward a value 
with a small temperature dependence. This is in agreement with our remark about the behavior of Inr as function 
of N and T in the paragraph above. We remark en passant that if following ^lOj and [7] we had ipwidth < "01 then 
A = l+CA^'''"'"''^''',with some constant C, and should accordingly converge towards one, with unfortunately extremely 
slowly decaying corrections. 

The quality of our data does not allow a trustable determination of the Kurtosis of the distribution G — 
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(Inr — Inr)"'/ (Inr — In r)-^ . For N = 256 and 512 the value of the Kurtosis is strongly affected by a few rare 

disorder samples with large values of Inr. Restriction the analysis to = 64 and 128, one obtains values of G « 5 
with little or no temperature dependence, but finite size effects. 

In conclusion, we have measured the probability distribution of the equilibrium relaxation time of the Sherrington- 
Kirkpatrick model with binary couplings, for a range of system size and temperature. We checked that our estimates 
are free of thermal noise. The data confirm that both the average and the width of the probability distribution of 
Inr scale as N-^^^ in the spin glass phase. As a last remark, we note that the long tail of the distribution of r may be 
a hidden source of severe problems in numerical simulations with the widespread practice of using the same number 
of Monte Carlo iterations for all disorder samples. This has been notices several time already, but is may be worth 
repeating. 
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